function [best_Zi, best_acc, Fmin, stage] = Bi_BPSO(num_simple, num_iter, num_feature, data_x, data_y, num_lable)
% num_simple:初始化粒子群样本数量;num_iter:搜索迭代次数;num_feature:每个样本包含的特征个数
% data_x代表聚类样本数据,data_y代表聚类样本的标签;num_lable代表聚类样本的标签的种类数
% best_Zi:准确率最高的样本的特征选择情况

% 初始化粒子群样本sample_x:num_simple*num_feature
X = rand(num_simple, num_feature);
% 将X转换成Z
r1 = ones(num_simple, num_feature) * 0.5;
Z = zeros(num_simple, num_feature);
Z(r1 < X) = 1;
Z(r1 >= X) = 0;
pbestx = X;
% 计算分类准确率F:num_simple*1
F = evaluate(Z, data_x, data_y, num_lable);
pbestf = F;
% 计算gbest,准确率最大的对应的索引
[fmin, gbest] = max(pbestf);
fprintf("iter 0 feature: pbestz(%d,:) fmin: %.4f\n",gbest, fmin);
num1 = 0;
Fmin = zeros(num_iter, 1);
stage = zeros(num_iter, 1);
for i = 1 : num_iter
    s = zeros(num_simple, num_feature);
    % pso更新下一步的位置,这里可以设置一下超过搜索范围的就设置为边界
    r = rand;
    if r > 0.5
        x = pbestx;
    else
        x = normrnd((pbestx + pbestx(gbest, :)) ./ 2,abs(pbestx - pbestx(gbest, :)));
    end
    r2 = rand(num_simple, num_feature);
    s(r2 < x) = 1;
    s(r2 >= x) = 0;
    presentf = evaluate(s, data_x, data_y, num_lable);
    % 更新每个单独个体最佳位置
    % pbestx(presentf > pbestf, :) = x(presentf > pbestf, :);
    % a reinforced memory strategy
    pbestx(presentf > pbestf, :) = 0.5 * (x(presentf > pbestf, :) + s(presentf > pbestf, :));
    pbestx((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)),:) = 0.5 * (x((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)), :) + ...
        s((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)), :));
    pbestf(presentf > pbestf, :) = presentf(presentf > pbestf, :);
    pbestf((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)),:) = presentf((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)), :);
    Z(presentf > pbestf, :) = s(presentf > pbestf, :);
    Z((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)),:) = s((presentf == pbestf)&(sum(Z, 2) > sum(s, 2)), :);
    % uniform combination (UC)
    pc = 0.2 / (1 + exp(5 - num1));
    for i1 = 1 : num_simple
        r3 = rand;
        if r3 < pc
            k = randi(num_simple, 1, 1);
            u_n = ceil(pc * num_feature);
            for i2 = 1 : u_n
                l = randi(num_feature, 1, 1);
                pbestx(i1, l) = pbestx(k, l) + rand;
            end
        end
    end
    % 更新所有个体最佳位置
    [fmin1, gbest] = max(pbestf);
    if fmin1 > fmin
        num1 = 0;
    else
        num1 = num1 + 1;
    end
    fmin = fmin1;
    fprintf("iter %d feature: pbestz(%d,:) fmin: %.4f\n",i, gbest, fmin);
    Fmin(i, 1) = fmin;
    stage(i, 1) = sum(Z(gbest, :));
end
best_Zi = Z(gbest, :);
best_acc = fmin;
end